Wave chaos as signature for depletion of a Bose-Einstein condensate 
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We study the expansion of repulsively interacting Bose-Einstein condensates (BECs) in shallow 
one-dimensional potentials. We show for these systems that the onset of wave chaos in the Gross- 
Pitaevskii equation (GPE), i.e. the onset of exponential separation in Hilbert space of two nearby 
condensate wave functions, can be used as indication for the onset of depletion of the BEC and the 
occupation of excited modes within a many-body description. Comparison between the multicon- 
figurational time-dependent Hartree for bosons (MCTDHB) method and the GPE reveals a close 
correspondence between the many-body effect of depletion and the mean-field effect of wave chaos 
for a wide range of single-particle external potentials. In the regime of wave chaos the GPE fails 
to account for the fine-scale quantum fluctuations because many-body effects beyond the validity 
of the GPE are non-negligible. Surprisingly, despite the failure of the GPE to account for the de- 
pletion, coarse grained expectation values of the single-particle density such as the overall width of 
the atomic cloud agree very well with the many-body simulations. The time dependent depletion of 
the condensate could be investigated experimentally, e.g., via decay of coherence of the expanding 
atomic cloud. 
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I. INTRODUCTION 



The workhorse for describing the non-equilibrium dy- 
namics of Bose-Einstein condensates (BECs) of ultracold 
gases is the Gross-Pitaevskii equation (GPE) (for a re- 
view see e.g. Ref. 1,2). Replacing the true many-body 
wave function by a single-particle orbital for the macro- 
scopically occupied condensate (particle number N 3> 1) 
results in an equation of motion that belongs to the class 
of nonlinear Schrodinger equations (NLSE). The GPE 
provides an appropriate starting point to investigate the 
underlying many-body system on the mean-field level. 
Effects beyond the GPE have been observed in BECs, for 
example in optical lattices with deep wells and small oc- 
cupation numbers per site.' 5 Other finite-number conden- 
sate effects include the demonstration of atom-number 
squeezing 4-6 and of Josephson junctions in a double 
well. Meanwhile, progress has been made in exploring 
the time-dependent many-bosons Schrodinger equation. 
One approach is the multiconfigurational time-dependent 
Hartree for bosons (MCTDHB) method which is a nu- 
merically efficient and, in principle, exact method for 
the time-dependent many-body problem. 10-12 In prac- 
tice, limitations are imposed by the finite yet large num- 
ber of configurations (millions) and orbitals (tens) that 
can be handled. 

We investigate repulsively interacting BECs after release 
into shallow one-dimensional (ID) potentials. The Bose 
gas is dilute and, initially, practically all particles are in 
one single-particle state. The external potential is weak 
compared to the single-particle energy. Comparison be- 



tween the MCTDHB method and the GPE for the expan- 
sion of the BEC provides detailed insights to what extend 
the GPE is capable of describing the condensate dynam- 
ics and may be capable of mimicking excitations out of 
the condensate state. One case in point is our recent 
observation of true (physical) wave chaos in the GPE, 1,! 
as opposed to numerical chaos 14 due to discretizations 
which has been exploited to study e.g. thermalization in 
the Bose-Hubbard system at the mean-field level. 15 Two 
wave functions nearby in Hilbert space arc exponentially 
separating from each other, as measured by the L 2 norm. 
Chaotic wave dynamics within the GPE is a mathemat- 
ical consequence of the non-integrability resulting from 
the interplay between the external (one-body) potential 
and the nonlinearity which replaces the inter-particle in- 
teractions. Its physical implications are, however, less 
clear as the original many-body Schrodinger equation is 
strictly linear and, thus, regular and non-chaotic. Previ- 
ously, a connection between chaotic dynamics within the 
GPE and growth in the number of non-condensed parti- 
cles has been made for time-dependent external driving 
which can be seen as a source of energy. 16 In the present 
study the external potentials are time independent such 
that the total energy is conserved. While wave chaos 
is likely associated with instabilities known for dynam- 
ics in periodic potentials (see e.g. Ref. 17) it is a much 
more general effect since it occurs in a broad variety of 
potentials ranging from harmonic oscillators with defects 
to periodic and disordered potentials. The aim of the 
present paper is to shed light on the physical meaning of 
wave chaos in the GPE for time evolution of BECs. By 
comparing the mean-field GPE and many-body MCT- 
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DHB dynamics we show that the build-up of random 
fluctuations within the GPE is directly related to the de- 
pletion of the condensate, i.e., the population of excited 
orbitals out of the condensate state. In this respect, the 
mean-field description surprisingly contains some infor- 
mation on many-body effects. Chaos on the mean-field 
level may open the door to identify depletion and ther- 
malization and the comparison with MCTDHB may shed 
light on how such a thermalization process takes place in 
a quantum many-body system. 

The outline of the paper is as follows. After first intro- 
ducing the system under investigation in Sec. II we briefly 
review the mean-field GPE and the many-body MCT- 
DHB method and identify relevant observables (Sec. III). 
The initial state whose dynamics we study upon release 
from the initial trapping is discussed in Sec. IV. We 
present numerical results for the dynamics in Sec. V fol- 
lowed by conclusions and remarks. 



II. SYSTEM UNDER INVESTIGATION 

We consider in the following a system of N bosons in- 
teracting via a pseudo-potential which captures the scat- 
tering dynamics of the real interaction potential in the 
limit of small wave numbers k — > 0. In a ID system with 
tight transverse harmonic confinement with oscillator fre- 
quency uj r the pseudo-potential is given by the contact 
interaction giT>8(x — x') with 
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(1) 



where a s is the 3D scattering length provided that a s <C 
y/H/mu) r such that the scattering can still be regarded 
as a 3D process. 18 The dynamics of the bosonic system 
is then determined by the many-body Hamiltonian (in 
second quantization) 



H = I dx — d x ^{x,t)d x '4){x,t) 



+ I dx V(x)ft(x,t)4>(x,t) 

dx if)* (x,t)^' (x,t) , 4>(x,t)'il)(x,t). 



9m 
2 



(2) 



The field operators fulfill the usual commutation rules for 
bosons. We study in the following the expansion of a Bose 
gas that is initially trapped also longitudinally (i.e. in the 
direction of expansion) by a harmonic potential with fre- 
quency ujq. These initial conditions serve to define char- 
acteristic sca les for l ength, time, and energy. We use the 
units Iq = y/h/mujQ for length, to = for time, and 

e = hujQ for energy. For a trap with ui r — lix x 70Hz and 
luq = 2tt x 5.4Hz used in a recent experiment on Ander- 
son localization 19 our units take on the numerical values 
Iq = 4.6/im and to = 29.47ms. 

Upon release from the trap, the particles move in an ex- 
ternal potential V(x) which we specify in the following 
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FIG. 1: (Color online) The initial (t < 0) harmonic trapping 
potential of radial frequency to r and longitudinal frequency 
uo (black line). At t — the longitudinal harmonic trap is 
switched off (the radial trap is preserved). Simultaneously, 
either a periodic (dashed red) or a disorder (blue solid) po- 
tential is switched on. The system expands for t > in the ID 
external potential. The length scale Iq corresponds to 4.6^m. 



to be a periodic potential of the form 



2tt 



V(x) — VaCos [ —x 



(3) 



with I = 0.54811?q (corresponding to I ~ 5.8£ with 
£ = h/y/Amfj, being the healing length and fi the chemical 
potential after release from the trap for N = 1.2 x 10 4 
87 Rb atoms) and varying potential amplitude Va- The 
periodic potential is realized in experiments by crossed 
laser beams in linear polarization along the same axis. 
For a realistic laser wave length tuned out of resonance 
with the 87 Rb 55 -» hP transition, \ L = 810nm, the 
above potential period of I corresponds to two linearly 
polarized crossed beams enclosing an angle of 8 0.1-7T. 
Alternatively, we also consider Gaussian correlated dis- 
order potentials Vd(x) of comparable strength with 
(Vd(x)) = 0, the variance Va = \J (V d 2 (x)), and the corre- 
lation length a. The interplay between the inter-particle 
interaction and the external potential plays a key role for 
chaotic dynamics resulting from non-integrability. 
We investigate in the following the dynamics of the ex- 
panding Bose gas in the mean-field approximation within 
the GPE and compare to the corresponding many-body 
dynamics within the MCTDHB method. 



III. METHODS 

A. Gross-Pitaevskii equation 

In the mean-field approximation the existence of a 
macroscopic occupation of one state is assumed such that 
the expectation value {tp(x,t)) = ip(x,t) takes on finite 
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values and can be treated as the classical field describ- 
ing the dynamics of the BEC. Further, requiring that the 
expectation value of the product of four field operators 
factorizes 



(^(x,t)^(x,t)^(x,t)^(x,t)} = |V0M)| 4 , (4) 

i.e., assuming second-order coherence, one arrives to- 
gether with 



ih 



&ip(x,t) 6(H) 



dt 



5ip*(x,t) 



(5) 



and H from Eq. 2 at the GPE 



ih 



dip(x, t) 
dt 



h 2 d 2 



2m dx 2 
gm\^(x,t)\ 2 '4)(x,t) 



ip(x,t)+V(x)ip(x,t) 



(6) 



Renormalization of the particle density to 
J dx\ip(x,t)\ 2 = 1 leads to the explicit dependence 
of the nonlinearity on the particle number N: 



ih 



chjj(x, t) 
dt 



h 2 d 2 



tp(x,t) + V(x)ip(x,t) 



2m dx 2 
g 1D N\xl>( X> t)\ 2 ij(x,t) 



(7) 



Since N enters as a parameter into the effective non- 
linearity, the GPE predicts the same dynamics for dif- 
ferent N as long as the product gmN is kept constant. 
In the limit iV — > oo with gmN — const, the (time- 
independent) GPE is expected to give exact results for 
the many-body system (at least for the ground state in 
three dimensions 20 ). 

The numerical simulations presented in the following are 
performed for N = 1.2 x 10 4 ultracold 87 Rb atoms in a 
cigar-shaped trap with frequency uj t (see Sec. II). The 
associated nonlinearity is 



go = 2huj r a s N ps 390eo^o- 



(8) 



Note that the rather high numerical value of go is due 
to the explicit inclusion of the number of particles N 
and does not contradict the assumption of weak inter- 
actions. Nevertheless, the interaction strength is suffi- 
ciently strong such that in the presence of an external 
potential depletion and fragmentation of the condensate 
may occur. 

To propagate the GPE, we use a finite element discrete 
variable representation (DVR) to treat the spatial dis- 
cretization (see e.g. Ref. 21,22). The propagation in time 
is performed by a second-order difference propagator (for 
details see Ref. 13 and references therein). 



B. Multi-configuration time-dependent Hartree for 
bosons (MCTDHB) method 



the condensate. Briefly, the many-body wave function is 
taken as a linear combination of time-dependent perma- 
nents 



{«} 



(9) 



where \ft\t) corresponds to states with occupation num- 
bers ft = (ni, njvf) and M is the number of single- 
particle orbitals. The sum runs over all sets of occu- 
pation numbers {ft} which fulfill ./V = X^i 71 *- ^ ne 
limit M oo the ansatz Eq. 9 gives the exact many- 
body wave function. MCTDHB efficiently exploits the 
fact that ultracold atoms may occupy only few orbitals 
above the condensate state. By dynamically changing the 
expansion amplitudes Cft(t) and the orbitals {$fc(r, t)}, 
even large many-body systems can be treated accurately. 
MCTDHB involves the solution of coupled linear differen- 
tial equations in C^if) and coupled nonlinear differential 
equations in {<6fe(r, t)}. The MCTDHB equations of mo- 
tion reduce in the case of M = 1 to the GPE Eq. 7 [with 
nonlinearity g\o{N — 1), the difference between TV and 
N — 1 can be neglected in the limit of large N] . 
Within the MCTDHB method kinetic operators are 
treated via a fast Fourier transform which is equiva- 
lent to an exponential DVR. 23 The nonlinear differen- 
tial equations in the orbitals are propagated via a 5 th or- 
der Runge-Kutta algorithm. The linear differential equa- 
tions for the amplitudes Cft{t) are propagated via a short- 
iterative Lanczos algorithm. The propagations are par- 
allelized using openMP and MPI. The MCTDHB orbital 
equations of motion are stiff differential equations for the 
interaction strength considered in this work. The appli- 
cability of the Runge-Kutta algorithm has been checked 
by comparing to an especially adapted integrator 24 ' 25 
ZVODE which relies on the Gear-type backwards differ- 
entiation formula for stiff ordinary differential equations 
and gives the same results as the faster Runge-Kutta al- 
gorithm. Numerical checks give a very good agreement 
between the initial and the backwards propagated den- 
sity per particle. The difference is of the order of 10 -4 
and less (in units of Iq 1 )- 



C. Observables 

The simplest and most important benchmark ob- 
servable for a comparison between the mean-field and 
the many-body dynamics is the single-particle density. 
Within the MCTDHB method the density is given by 



p(x,t) = N j dx 2 ,...,dxN ^*{x,x 2 ,---,XN;t) 
x y(x,x 2 ,-~,x N ;t) 

M 

= J2 Pm,n( t Wm(x,-t)^n(x,t), (10) 



The MCTHDB method 10 ' 11 allows one to describe 
many-body effects beyond the mean-field description for 



where the elements p mn {t) are readily accessible as a 
combination of the amplitudes C^(t) and the correspond- 
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ing occupation numbers contained in ft. Upon diagonal- 
ization of Eq. 10 the natural orbitals &f°(x, t) and their 
occupation numbers nf°(t) are obtained as: 

M 
i=l 

In the presence of a BEC the occupation of one state has 
to be "macroscopic" 26 (of order N) . In the following we 
denote this condensate state as $fS\(x,t) and its occu- 
pation as nfS[(t). All other states $f°(x,t) with i > 1 
are referred to as excited states. The Fourier spectrum 

M 

f>(k,t)= Pm,n(*)$m(M)$n(M), (12) 

m,n— 1 

is obtained by Fourier transforming the orbitals $ n (x,t) 
to give $> n (k,t). In the limit of a long-time expansion 
of the BEC in free space when the initial interaction en- 
ergy is converted into kinetic energy p(k, t) corresponds 
to the experimentally observed momentum distribution. 
Within the GPE, p(k, t) is given by the absolute square 
of the Fourier transform of the condensate wave function 

m,t)\\ 

We utilize coherence as measured by the normalized two- 
particle correlation function 2 ' 28 

(2 ), , , .x _ P {2) (x' 1 ,x' 27 x u x 2 ;t) 

g y >(x 1 ,x 2> Xi,X2 ) t) — = 
V p(xi,t)p(x 2 ,t)p(x' 1 ,t)p(x' 2 ,t) 

(13) 

to analyze the breakdown of the GPE on the length scales 
of the random fluctuations which develop in the wave 
function in the regime of wave chaos. In g^ the reduced 
two-body density matrix 

p ( - 2 \x' 1 ,x' 2 ,x u x 2 ;t) = N(N - 1) 

x J dx 3 , ...,dx N ^*(x' l7 x 2 ,x 3 , ...,x N ;t) 

x^(xi,x 2 ,X3,...,x N ;t) (14) 

enters. For a fully second-order coherent system g^ 
fulfills \g^(x' 1 ,x 2 ,x 1 ,x 2 - 1 t)\ = 1. Within the GPE the 
reduced two-body density matrix is a product of one- 
body wave functions. Thus, |g( 2 )| = 1 for all times, 
i.e., full second-order coherence is a generic feature of 
the GPE. In the many-body case the departure of If/ 2 -*! 
from unity gives a measure for how well the system is 
described by a single-orbital product state and how cor- 
related d^ 2 '! > 1) or anticorrelated (Iff' 2 -'! < 1) the mea- 
surement of two coordinates is. (Anti-)Correlation indi- 
cates the degree of fragmentation in the system. 
As a measure for wave chaos, i.e. the build-up of random 
local fluctuations on the length scale comparable to that 
of the external potential, we have introduced 13 the Lya- 
punov exponent characterizing the exponential increase 
of the distance in Hilbert space of two initially nearby 
GPE wave functions tpi j2 (x, t). The distance is measured 



by the L 2 norm 

d {2 \t) = \ J dx |V>i(2,t)-V>2(z,*)| 2 

= 1 -Re (J dx ipl{x,t)i>2{x,t)\ . (15) 

The distance function takes on values d^ € [0, 2] and 
is 1 for orthogonal wave functions. In terms of d^ 2 ', the 
Lyapunov exponent indicating chaos is given by 

A = I lim lim - In ( ■ (16) 

2 t-foodw (o)^o t \d( 2 H0)J V ; 

e?( 2 ) is invariant for unitary time propagation of lin- 
ear systems: if ipx(x,t) and ip 2 (x,t) would be solu- 
tions of the linear Schrodinger equation, d^ would be 
constant. Similarly, inserting instead of tpi(x,t) and 
ip 2 (x,i) two many-body wave functions ^(xi, Xjy, t) 
and ^(xi, ■■■,XN,t), and integrating over all spatial co- 
ordinates gives a constant of motion. By contrast, con- 
struction of a reduced one-particle wave function from an 
initial A-body state of system 1 by (see e.g. Ref. 29) 

^OM) = (*i(JV-i)$(M)|*i(i\0) 

= J dx 2 ...dx N *i(x 2 , ...,x N ,t) 

x^i(x,x 2 ,...,x N ,t), (17) 

where |^i(A)) denotes a many-body state with N parti- 
cles leads to a many-body measure analog to the d^ (t) 
function that is not conserved as a function of time. This 
can be seen by inserting Eq. 17 into Eq. 15, taking the 
time derivative of d^ 2 \ and using the many-body Hamil- 
tonian. In the time derivative of d^ 1 contributions orig- 
inating from the kinetic energy and the external poten- 
tial V(x) cancel, while contributions from the interaction 
term lead to d[d^ 2 '{t)]/dt 7^ 0. The tracing out of unob- 
served degrees of freedom leads to the violation of the 
distance conserving evolution. In the particular case of 
the GPE the nonlinearity present can cause exponential 
divergence of d^ 2 \ 

It is now of key importance to relate the behavior of 
the d^ function within the GPE to properties of the 
unitary time evolution of the underlying many-body sys- 
tem, as described by the MCTDHB method. The work- 
ing hypothesis is that the random fluctuations devel- 
oping within the GPE are the signature for its failure 
to properly account for the depletion of the condensate 
state, i.e. excitation of the BEC during expansion in 
an external potential. In turn, within the MCTDHB 
method the population of all natural orbitals beyond 
that describing the condensate should grow. While at 
t = the MCTDHB method and the GPE closely agree 
which each other with only one natural orbital occupied, 
f^°(0) = <°(0)/A « 1 (see next Sec. IV), with in- 
creasing time all other occupation numbers nf (i > 1) 
should increase. In the following we study the dynamics 
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of N = 10 3 to 10 5 particles for which the ground state 
densities closely agree with each other (see Fig. 2). For 
N = 10 4 and TV = 10 5 only M = 2 orbitals allow a nu- 
merically feasible number of configurations (N c = 10 4 + 1 
to N c — 10 5 + 1, respectively). Already adding one 
more orbital (M = 3) leads to a configuration size of 
N c = 50,015,00f for N = 10 4 which may be at the 
border of feasibility and requires a massive paralleliza- 
tion over a large number of processors. The system with 
N = 1 5 and M = 3 resulting in N c « 5 x 10 9 is out of 
reach for the current implementation of the MCTDHB 
method. For N = 10 3 a number of orbitals up to M = 3 
is numerically feasible and allows to quantify the effect 
of adding one more orbital to the case M = 2. Due to 
the numerical limitations we focus on the early stages of 
the depletion process when the depletion is still relatively 
weak nf°(t) > 0.95. 

As a quantitative measure for the depletion process we 
introduce the state entropy for a general many-body state 



S N (t) 



$>f (i)lnnf°(;) 



(18) 



where 



:NO 



(*) 



,NO/ 



(t)/N and nf (t) are the occupa- 



tion numbers of natural orbitals obtained by diagonal- 
izing the reduced one-body density at each time (see 
Eq. 11). For the initial conditions used in the present 
study we have Sjv(t) > Sn(0) « 0. Sjv(i) measures 
the entropy of the occupation numbers of the many-body 
system along distinct quantum states. Within the GPE, 
where nf °(t) = 1 and nf-^\(t) = 0, SW(i) remains strictly 
zero. Deviations of Spf(t) from zero within the many- 
body theory thus mark the breakdown of the mean-field 
theory within the GPE. In the following we will focus 
on the time evolution of Eq. 18 and investigate the time 
scale of depletion, i^, defined by the occurrence of an 
abrupt change of Sn from Sn ~ to Sn > 0. We 
associate this observable with the onset of exponential 
growth, t e , of S 2 )(t) within the GPE. t e is determined 
from the crossing-over point between the free-space ex- 
pansion behavior of S- 2 \t) to the exponential increase. 
td is implicitly dependent on N through the degree of co- 
herence of the condensate. The larger N, the smaller the 
depletion (nf°, i > 2) at the same time. Consequently, 
the depletion time is size dependent td(N) (see below). 



IV. THE INITIAL STATE 

The initial state of the bosonic gas corresponds to the 
ground state of the harmonic trap. For this ground state 
the GPE with nonlinearity go predicts a BEC in the 
Thomas-Fermi regime. Applying both the MCTDHB 
method and the GPE to the same system requires a care- 
ful choice of system parameters, in particular the parti- 
cle number N. While the validity of the GPE calls for 
the limit of large N — > oo, such a case is numerically 
prohibitive for the MCTDHB expansion Eq. 9. Since 
the GPE results are invariant for varying N but fixed 



.9o = gir>N we adjust the particle number such as to 
remain in the Thomas-Fermi limit of the longitudinally 
trapped BEC (see Fig. 2). In such a way it is assured 
that discrepancies between the GPE and the MCTDHB 
method during the time evolution are not caused by in- 
compatible initial conditions. 

The ground state of an interacting system of bosons 
trapped by a harmonic potential is governed by three 
length scales: the characteristic length of the har- 
monic trap Iq, the mean inter-particle distance r s = 
nT x with n the particle number per unit length, and 



lx = a measure for the zero-point fluctuations (or 

anti-correlation length) of the repulsive two-body delta- 
function interactions of strength <?id- The regimes ob- 
tained range from a non-interacting "Gaussian" shaped 
BEC, over a Thomas- Fermi BEC, to a strongly interact- 
ing fermionized Tonks- Girardeau gas. 30 The presence or 
absence of a BEC is determined by the ratio 



7 : 



(19) 



referred to as the Lieb-Lininger parameter. 31 If Is is much 
larger than the inter-particle spacing r s the particles are 
allowed to occupy the same state and form a BEC. The 
condition for the presence of a BEC thus is: 



7< 1. 



(20) 



In order to distinguish between a Gaussian and a 
Thomas-Fermi BEC the harmonic oscillator length Iq 
must be included. The regimes are controlled by the 
parameter 30 



lo 



a = — = 



gig 



(21) 



In our case a <C 1 (a = go/N with go from Eq. 8), 
i.e. Iq <C Is- In this case the system is in the Thomas- 
Fermi regime if in addition Iq 3> r s which is equivalent 
to the statement that the Thomas-Fermi length Ltf = 
-y/2/i/mojQ corresponding to the width of the BEC in 
the Thomas- Fermi regime is much larger than l . The 
condition 3> r s implies N 3> a" 1 for the Thomas- 
Fermi limit to hold. 30 For all systems with N > 1000 in 
Tab. I the criteria N 3> or 1 and 7 <IC 1 are well fulfilled 
and, indeed, the many-body ground state density takes 
on the Thomas- Fermi shape (see Fig. 2). The density 
is practically indistinguishable from the GPE prediction. 
For comparison, we also show in Fig. 2 a system with 
N = 100 for which the criterion of a Thomas-Fermi BEC 
is only marginally fulfilled because 7 sa 1 and deviations 
become apparent. 
The requirement of large N places a severe limit on the 
number of orbitals that allow for a numerically feasible 
configuration space. Convergence in the orbital number 
is controlled by the occupancy n^f of the least occupied 
state. While for the ground state calculations is 
sufficiently low, we expect this number to rapidly increase 
during expansion since strong depletion may occur which 
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FIG. 2: (Color online) The initial state of the BEC in the 
harmonic trap with ujq — 2n x 5.4Hz and Iq w 4.6/^m. Black 
solid line: GPE. Red dotted line: MCTDHB with TV = 1000, 
M = 3. Blue dashed line: MCTDHB with TV = 100, M = 3. 
The two local maxima for TV = 100 are due to depletion of 
the condensate in the initial state and indicate deviations from 
the Thomas-Fermi limit. 

will spread over a large number of orbitals. 
We, therefore, expect only the onset of depletion to be 
quantitatively reliable while the occupation numbers of 
excited orbitals can be considered to be an indication 
of the excitation process as the orbital expansion ceases 
to converge (M > 3 time-dependent orbitals would be 
needed) with increasing propagation time. 



V. NUMERICAL RESULTS 

We first consider the expansion of the BEC initially 
formed within the harmonic trap (Fig. 2) which is, after 
the release, subject to a periodic potential (Eq. 3) with 
I = 0.54811^o and Va = 0.2e with e the total energy per 
particle. After the release an explosion-like process takes 
place: the interaction energy is rapidly transformed into 
kinetic energy. In free space the cloud expands keeping 
its Thomas-Fermi shape with the characteristic length 
increasing in time. 32 This process is modified by the 
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0.39 
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10 4 
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0.039 


0.9997 


0.33 x 10~ 3 


10 5 
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0.0039 


0.99997 


0.34 x 10" 4 



TABLE I: Parameters of the many-body systems trapped in 
the harmonic oscillator at t = (Fig. 2): Particle number 
TV, number of orbitals M. The interaction strengths giu cor- 
respond to constant go = giuN (Eq. 8). For the definition 
of a see Eq. 21. Highest occupation number n^°, smallest 
occupation number nff. The parameter 7 (Eq. 19) fulfills 
7< 1. 
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FIG. 3: (Color online) (a) The density of an expanding cold 
atom cloud in a periodic potential within the MCDTHB for 
TV — 10 5 at t = 3io- The inverted parabola in blue dashed 
line indicates the shape of an expanding cloud in the absence 
of the potential. The parameters of the periodic potential 
Eq. 3 are / = 0.54811/o and Va = 0.2e, where e is the energy 
per particle on the mean-field level, (b) The corresponding 
Fourier spectrum again compared to its form for free expan- 
sion (blue dashed parabolic curve). The vertical lines corre- 
spond to ±k L = ±4.434/o \ 



presence of the periodic potential. Practically immedi- 
ately the density is modulated by standing waves with 
the same spatial periodicity as the potential. The local 
maxima of the density coincide with the local minima of 
the potential and lead to an increase of kinetic and in- 
teraction energy at cost of potential energy. 
As soon as the Fourier spectrum is sufficiently broad, 
inelastic processes set in. At that point a fraction of 
particles have acquired a momentum k > k^ = 
with vl the Landau velocity which leads to excitations 
of phonons, i.e. friction of superfluid flow. For a homoge- 
neous system vl is given by vl = U with n the par- 
ticle density and <?id the interaction strength (prefactor 
of contact interaction). An estimate for the Landau ve- 
locity in case of an inhomogeneous system can be derived 
from the chemical potential fj, which contains an average 
over the small scale fluctuations of the density. For a 
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homogeneous system /j = ngm- We thus use vl = \ 
to determine the numerical value of k^. At f ~ 3io the 
width Afc of the Fourier spectrum is approximately as 
large as kh and we observe the development of strong 
density modulations [Fig. 3 (a)]. These spatial density 
modulations go hand in hand with reduced density in the 
Fourier spectrum near ±k^ since those particles lose their 
momentum by phonon excitations [Fig. 3 (b)]. Friction 
leads to the separation of a strongly fluctuating central 
part of the density from its fast tails [Fig. 3 (a)]. The 
tails expand nearly freely and are modulated by the po- 
tential. We point out that this process is fully accounted 
for within the GPE (i.e., the system remains condensed) 
since it gives practically the same density and spectrum 
for t = 3t as MCTDHB in Fig. 3. 

For longer times we have previously observed for this sys- 
tem signatures of wave chaos: 13 two nearby effective one- 
body wave functions ipi(x,t) and ifj2(x,t) (with initially 
large overlap) propagated by the GPE become orthogo- 
nal to each other after an exponential increase in distance 
in Hilbert space [see Fig. 4 (c)]. This orthogonality re- 
sults from the build-up of random local fluctuations in 
the wave functions on length scales comparable to the 
period of the potential. It is now instructive to compare 
the growth in g?( 2 ' within the mean-field description with 
the growth of S N {t) (or depletion) within the MCTDHB 
method which the GPE cannot follow. For vanishing 
potential V(x) = we find that the explosion-like expan- 
sion with a rapid transformation of interaction energy to 
kinetic energy does not lead to depletion of the conden- 
sate [Fig. 4 (b) dashed line]. The GPE accounts for the 
expansion dynamics since Sjsr(t) remains approximately 
zero as a function of time. For vanishing potential the 
GPE is integrable 33 such that saturates after a short 
universal increase [Fig. 4 (c) dashed line]. For periodic 
potentials (Fig. 4) we find over a wide range of poten- 
tial strengths (0.04e < Va < 0.2e) that the exponential 
separation on the mean-field level and the depletion on 
the many-body level correlate well with each other. To 
extract the depletion time td and the rate of depletion r\ 
we fit the curves to functions of the form 
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FIG. 4: (Color online) (a) The potential V(x) in units of the 
energy per particle e. The amplitude Va is varied succes- 
sively, the period I is kept constant I = 0.54811/n- (b) The 
state entropy S N (t) within MCTDHB and (c) d (2) function 
of two close wave functions within the GPE expanding in the 
periodic potentials of (a) . The parameter Va is indicated next 
to the curves in units of the energy per particle e. The many- 
body system corresponds to a particle number of N = 10 4 
and an orbital number of M = 2. The onset of exponential 
growth t e as well as the onset of depletion td are marked by 
arrows, td corresponds to a time of « 60ms. Note that Sn is 
strictly zero within the GPE. 



with fit parameters c, a, and td (in Fig. 4 td as well as t e 
are marked for Va — 0.2e). in the above equation is 
the Heaviside step function. We introduce the depletion 
rate 77 as 



77 = 



(23) 



The depletion rate 77 is equal to the slope of Sj^(t) 
t = td + to, i.e. after the abrupt increase of S^(t) at td- 
77 has the same dimension as the Lyapunov exponent A. 
An increasing Lyapunov exponent A with increasing Va 
goes hand in hand with an increasing 77 (Fig. 4). Up to 
a constant numerical factor (sa 10) 77 closely follows A 



at 



as a function of Va and even a threshold-like behavior is 
present (Fig. 5, compare also to Ref. 13 Fig. 7). We note 
that both td and 77 are sensitive to the specifics of the fit 
function Eq. 22 which results in an uncertainty of the fit. 
Since we are interested in the qualitative behavior here, 
the error bars in td and 77 are of no concern. 
The variation of td from the onset of exponential diver- 
gence of S 2 ^ in the GPE is comparable to the variation of 
td within the MCTHB method for different particle num- 
bers N [compare Fig. 4 (b), (c) with Fig. 6]. td increases 
with increasing N and approaches rj d = 6.48in which we 
obtain from the GPE by associating the onset of chaos t e 




V A /e 

FIG. 5: The depletion rate r\ as a function of Va as compared 
to the Lyapunov exponent A. rj and A are calculated from 
Sn (i) and d 1 - 2 ' , respectively, in Fig. 4. 

with the onset of depletion in the limit N — > oo. For rela- 
tively small N (N = 10 3 ) the MCTDHB simulations are 
also feasible for M = 3. Comparing M = 2 and M = 3, 
the threshold for depletion is only weakly dependent on 
the number of orbitals included: we obtain almost the 
same td for M — 2 and M = 3. 

Another important example is propagation in a disorder 
potential. We compare the results for periodic poten- 
tials to propagation in a disorder potential with parame- 
ters for which previously Anderson localization has been 
observed. 19 Within the GPE for Va = 0.2e and corre- 
lation length a = 0.7£ we obtain on a time interval of 
t G [6, 15] to an exponential increase with practically the 
same Lyapunov exponent as for the periodic potential 
with Va = 0.2e and I = 0.54811/ (not shown). Only 
initially the slope is reduced [Fig. 6 (c)]. Accordingly, 
we find within the many-body system that the curves for 
Sn (t) follow each other with a reduced rate for the disor- 
dered system [Fig. 6 (b)]. We point out that our results 
indicate a destruction of the BEC as signified by the oc- 
cupation of excited states during expansion in disorder 
potentials in the regime of Anderson localization. 
The onset of depletion of the condensate is mirrored in 
the fine scale oscillations of the density (Fig. 7) . Substan- 
tial deviations of the density within the GPE from the 
density obtained within the MCTDHB method emerge 
at different instants of time for different particle num- 
bers. For N — 10 4 deviations in the local fluctuations 
of the density emerge at t ~ 4to (see Fig. 7) monitored 
by Sjv(i) > 0. The occupation numbers are « 0.96 
and ~ 0.04. While the condensed part [given by 
n^°(t)\<f>%°(x,t)\ 2 } still closely follows the GPE predic- 
tion \ip(x, t)\ 2 , the total density p(x) shows smoothing 
of the local fluctuations near the center. This smooth- 
ing is due to excited atoms whose density partially fills 
in the local minima. For the system with N = 10 5 the 
picture is very similar except that the occupation of the 
excited state is lower at t = 4t , n^ ~ 0.003 instead 
of « 0.04. Spatially strongly localized structures 
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FIG. 6: (Color online) (a) The periodic and disordered po- 
tential with amplitude Va = 0.2e. The period of the periodic 
potential is / = 0.54811io ~ 5.8£. The correlation length of 
the disordered potential is a — 0.7£. (b) Onset of depletion 
within MCTDHB for different particle and orbital numbers 
for propagation in the periodic and disordered potential of 
(a). The numbers next to the curves indicate the particle 
number N. (c) Onset of chaos within the GPE. In both (b) 
and (c) the thick line corresponds to the disordered poten- 
tial. The onset of depletion td and exponential divergence i e 
is marked by arrows. 



emerge which the GPE overestimates as the transition 
into excited orbitals is forbidden. The initially spatially 
localized excitations spread over the entire system with 
increasing time. One can expect the fine scale structure 
of the density of the full many-body system to strongly 
differ from the prediction of the GPE. 
The discrepancies in the particle density go hand in hand 
with the breakdown of coherence as measured by the nor- 
malized two-particle correlation function (Fig. 8): 
directly quantifies the degree of second-order coherence 
in the system 27 ' 28 and furthermore can be measured in 
the experiment (see e.g. Ref. 34). In the regions of high 
density of excited atoms (near the local maxima of the 
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FIG. 7: (Color online) Particle density at t = 4to within 
MCTDHB for N = 10 4 and GPE. The condensate density is 
given by (t)\$^ (x, t)\' 2 , the density of excited atoms is 
determined by n§°(t)\<f>%°(x, t)\ 2 . 
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FIG. 8: (Color online) Two-particle normalized correlation 
function g^- 2 ' {x\,X2,xx,X2 1 t = 4to) for iV = 10 4 bosons as 
well as the density of the second natural orbital \cf>2 (x,t = 
4t )| 2 (inset) obtained from MCTDHB. White corresponds 
to no correlations and full coherence, red to correlation and 
blue to anti-correlation, i.e. reduced coherence. Note that 
g( 2 \xi, X2, Xi, X2\ t) is a real function (see Eq. 14) and is equal 
to unity within the GPE. 



second natural orbital) the two-particle coherence is lost, 
i.e. Q strongly differs from 1. The deviation of from 
unity indicates that the many-body state is no longer 
representable by a product of a single complex-valued 
function, i.e. the GPE ceases to be a valid description. 
This is a fingerprint of the emerging fragmentation of 



the many-body system. Roughly speaking, the pattern 
of g^(xi,X2,Xi,X2]t) is similar to the functional form 
of \4>2°(xi,t) x (j)2°( x 2, t)\- For longer time intervals our 
MCTDHB calculations indicate a destruction of the con- 
densate or at least a strong fragmentation. At t — 20io, 
e.g, the occupation of both orbitals is approximately 50% 
indicating that many more orbitals would be required for 
convergence. Nevertheless, current experiments indicate 
remarkable agreement with the prediction of the GPE 
for large scale observables such as the width of the atom 
cloud or the average position (see e.g. Ref. 19,35-37). For 
the width 



Ax = y/(x 2 ) - <x) 5 



(24) 



where (x n ) — J dxp(x)x n , we previously observed that 
within the GPE it is independent of wave chaos: two 
close wave functions ipi(x, t) and ^{x, t) develop random 
local fluctuations and become orthogonal. The width, 
however, for both tpi(x,t) and ip2(x,t) agrees. If we now 
compare the prediction for the width within the GPE 
and within the MCTDHB method, we observe the same 
trend. Despite its failure to account for the state entropy 
(Fig. 4 and Fig. 6) and the coherence properties (Fig. 8) , 
the GPE remains predictive in describing the expansion 
of a BEC in external potentials on longer time scales 
for coarse-grained observables such as the width and the 
average position, long after the random fluctuations pre- 
vent the prediction of fine scale structures in p(x). Up 
to now local small-scale fluctuations have not been inves- 
tigated experimentally because of the difficulty of (sub) 
fim resolution. While the fine scale structures of the wave 
function within MCTDHB have not fully converged for 
the small number of orbitals (M < 3 included in the 
simulation), the coarse-grained distribution remains es- 
sentially unchanged compared to the GPE [Fig. 9 (a)] 
and for varying N and M. We thus expect that the time 
dependence of the width of the full many-body system is 
well accounted for by the curves in Fig. 9 (a). The same 
excellent agreement we observe for the average over mo- 
menta k 2 as accessible in time-of-flight experiments [see 
Fig. 9 (b)]. The average over k 2 is determined via 



(k 2 



dk k 2 p(k,t) 



(25) 



and is proportional to the kinetic energy per particle. 
The GPE thus reproduces the mean kinetic energy of a 
highly excited system despite its failure to account for 
breakdown of coherence, fragmentation, and small scale 
fluctuations. The latter observation indicates that ther- 
malization may be within the realm of the GPE despite 
its failure to account for two-body scattering which is key 
to any thermalization process. 



VI. CONCLUSIONS 

By comparative simulations invoking the Gross- 
Pitaevskii equation (GPE) and the multiconfigurational 
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FIG. 9: (Color online) (a) The width A function of time 

as predicted by the GPE (solid black line) and the MCTDHB 
for N = 10 3 (red squares), N = 10 4 (green circles), N = 10 5 
(blue triangles), (b) Average of the square of momentum (k 2 ) 
(or mean kinetic energy) as a function of time, symbols as in 
(a). 



time dependent Hartree for bosons (MCTDHB) method 
we have uncovered a direct connection between wave 
chaos in the GPE and depletion of the occupation of 
a BEC during expansion in the presence of weak ex- 



ternal ID potentials. We have checked that this con- 
nection holds for a large class of external potentials in- 
cluding a harmonic potential with short-ranged pertur- 
bation, an aperiodic potential with incommensurate fre- 
quencies, and disordered and periodic potentials explic- 
itly discussed in this paper. This connection has far- 
reaching consequences: while the depletion and fragmen- 
tation process is an intrinsic many-body effect outside 
the realm of the GPE, the mean-field theory allows one 
to monitor its onset through the development of random 
local fluctuations. The measure for the random local fluc- 
tuations, S- 2 \ can be used to delimit the applicability of 
the GPE to approximate the many-body dynamics. On 
the many-body level the depletion process manifests itself 
through the loss of coherence as measured by deviations 
of g( 2 ) from unity. The connection between wave chaos 
and depletion indicates that a mean-field description may 
allow to identify and describe relaxation and thermal- 
ization processes (see e.g. Ref. 38-41), when the system 
reaches equilibrium for coarse-grained quantities while 
non- universal fluctuations about the equilibrium persist. 
The depletion process, the onset of which we have inves- 
tigated, can be experimentally studied provided a suffi- 
cient spatial resolution is achieved. Observables include 
higher-order coherence, i.e. deviations of from unity 
as measured e.g. in Ref. 34. It would be of consider- 
able interest to verify experimentally our predictions by 
exploring the fine-scale fluctuations and coherence prop- 
erties of expanding BECs in external potentials and thus 
gain deeper insight into the involved many-body effects. 
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